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1. Introduction 

When doing error analysis for numerical quadrature, achieving good uniform 
bounds on higher order derivatives of the integrand is paramount. As undergrad- 
uates become increasingly adept with programmable calculators, numerical inte- 
gration schemes such as Simpson's Rule and the Trapezoidal Rule take on a new 
relevance. Although it may be a rare calculus class that dwells overmuch on er- 
ror bounds for such schemes, this may be due as much to the perceived paucity 
of interesting examples for which decent error bounds are readily achievable as to 
the general weakness in algebraic skills necessary for the requisite understanding 
of inequalities. The purpose of this article, therefore, is to offer some interesting, 
non-trivial examples for which the error analysis, if not elegant, is at least simple 
enough to carry out in the classroom. 

2. The Sine, Cosine and Exponential Integrals 
The sine integral ([T| pp. 231-232] or pH pp. 503-504]) is given by 

Si(x) = f — dt. (2.1) 

The integrand of (|2.ip is a scaled version of the sine function, as it is referred to 
by those who work on problems involved with signal processing and reconstruction. 
Thus, sinc(f) = sin(irt) / (nt) for t ^ 0, and it is convenient to define sinO/0 = 1, 
thereby removing the removable discontinuity. For Simpson's rule, we require a 
good uniform upper bound on the absolute value of the fourth derivative of the 
integrand ([H p. 886] or [3 pp. 57-58]). Repeated application of the product rule 
for differentiation gives 

/ d\ 4 sint sini A cost 12 sin t 24 cost 24 sin t 

\dt) ~T~ = ~ + ~ T 3 ~ + ~ ¥~' 

which in this form is difficult to estimate, not only because of the complicated 
nature of the formula, but also because of the apparent difficulties when t is near 
zero. However, 

sini f 1 , v , 

= / cos(st) ds, (2.2) 

t Jo 
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and hence by Leibniz's rule for differentiating under the integral sign, 



sini 



s 4 cos(st) ds. 



Since | cos(s£)| < 1 for all real s and t, it follows that 



sini 



< 



s 4 ds = 



The inequality is sharp, as can be seen by substituting t = in 
manner, one can show that for all nonnegative integers k, 



s'mt 



< 



k ds = 



1 



1 



(2.3) 

(2.4) 
In a similar 

(2.5) 



with equality at t = when k is even. Instead of starting with (|2.2p . an analysis 
based on the equivalent representation 

can likewise be given, although in freshman/sophomore calculus classes it is less 
likely that students will be comfortable with the complex exponential. 
For the cosine integrals {[1, pp. 231-232] or pH pp. 503-504]) 



Cin(x) = 



1 — cos t 

t 



dt, 



Ci(x) = 7 + log a; + 



x cos t - 1 

t 



dt, 



we define the integrand to be zero when t = so that for all real t, 



and hence 



1 — cos t 
t 



1 — cos t 

t 



sin(st) ds, 



s 4 sin(st) ds. 



(2.7) 



(2. 



Since j sin(st)| < 1 for all real s and t, it follows that 



1 — cos t 



t 



< 



and in general, for all nonnegative integers fc, 



1 — cos t 



t 



< 



s ds = — , 
5 



s fe ds 



1 



(2.9) 



(2.10) 



with equality at t = when k is odd. 

The hyperbolic sine and cosine integrals ([U p. 231] or [7, p. 936]) 

rX sinht 



Shi(a;) 
Cinh(a;) 
Chi(a;) 



t 



■ dt, 



x 1 - coshi 
t 



7 + log x + 



dt, 

x cosh t — 1 



t 



dt, 



x>0 
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can be treated in much the same way. For example, for Shi one replaces sin by sinh 
and cos by cosh in (I2.1|) - (I2.3|) . The bound in (|2.5p is modified by the presence of 
an extra factor of coshi or | sinht according to whether k is even or odd. In either 
case, equality obtains when t — 0. 

The exponential integral ([H pp. 228-231] or Qjfl pp. 504-505]) 



E 1 (x) = J 



oo e - xt 



dt, x > (2.11) 

is another natural choice to study using this technique. In this case, it is easier to 
deal with the complementary expression 

f x 1 - e"* 

Ein(x) = / — dt, (2.12) 

Jo t 

since the latter defines an entire function. In view of the relationship ([TJ p. 228] 
or p. 40]) 

Ein(x) = logx + 7 + Ei(x), x > 0, 
there is no essential difference between the two. If we define the integrand of (|2.12[) 
to be zero when t = 0, then 

l ^r 1 'll'-' da - {i)" Lz f 1 = Jj- s) '' e -'<' s - (2 - 13) 

and so for t > and k a nonnegative integer, 

*i pi 



d \ k 1 - e"* 



s k e- sf ds< I sfcrfs =£-|-p ( 2 - 14 ) 



Jt J t 
with equality again when t = 0. 

3. How Good Are These Estimates in Practice? 

The estimates (|2.4j) and (|2.14p of the previous section are best possible in the 
sense that equality holds in each when t = 0. On the other hand, the corresponding 
derivatives (|2.3[) and (|2.13[) each tend to zero as t grows without bound, so it is clear 
that a uniform numerical bound for the entire range of t values is less than ideal. 
Nevertheless, it is instructive to see just how well our estimates hold up in practice. 
We confine ourselves here to a single example, the cosine integral, in which we test 
the inequality (|2.9p used to give an error estimate for Simpson's Rule against the 
error arising from an actual computation. 

Let / : R -» R be given by /(0) = and f(t) = (1 - cost)/* for t ^ 0. Then the 
cosine integral (|2.7p is given by 

Cin(x) = / f(t)dt. (3.1) 



For even positive integers n and real x > 0, define 

, n/2 > , . v n/2 — 1 > v ^ 

the approximation to the integral (|3.ip obtained by applying Simpson's Rule with 
n subdivisions of the interval [0,x]. The error is given by [51 p. 58] 



E n {x) := S n (x) - Cin(x) 



x 



/ (4) (?) 



180n 4 
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where £ = £(x) £ [0,x]. In view of (|2.9|) . we have the inequality 

\E n (x)\ <B n (x), where ^):=^. (3.2) 
To compare the estimate £?„ with the actual error E n , we let 

r ( x ) ■= Bn ^ = Bn ^ 

E n (x) S n (x) — Cm{x) 

be the ratio of the two quantities. Then R n is defined except for those points at 
which the actual error vanishes, and elsewhere \R n \ > 1. 

Using Maple's built-in cosine integral, it is possible, in principle, to compute E n 
(and hence R n ) to any desired precision. Table [1] gives an indication of the range 
of values of R n (x) for < x < 10. Tabular values were computed using Maple V 
Release 5 with a conservative working precision of twenty decimal places. Entries 
were then rounded to two decimal places to make comparison easier. Thus, for 
example, the error bound given by (|3.2[) arising from the estimate (|2.9j) exceeds the 
actual error by a factor of almost 3 when x — 5 and Simpson's Rule is used with 
n = 10, 100, or 1000 subdivisions. When x — 7, the corresponding overestimate is 
by a factor of less than 12. 



x \ n 


10 


100 


1000 


1.0 


2.55 


2.56 


2.56 


2.0 


1.54 


1.55 


1.55 


3.0 


1.43 


1.44 


1.44 


4.0 


1.75 


1.77 


1.77 


5.0 


2.77 


2.82 


2.82 


6.0 


5.63 


5.75 


5.75 


7.0 


11.31 


11.58 


11.58 


8.0 


10.66 


11.18 


11.18 


9.0 


6.97 


7.51 


7.51 


10.0 


5.59 


6.14 


6.15 



Table 1 . Values of R n (x) 



Figure Q] is a graph of Riq(x) for < x < 34, and makes it clear that there is 
nothing special about our choice of integer values for x in Table Q] We see that 
for < x < 20, 1 < Rio(x) < 15, i.e. the bound (|3.2[) is never farther than a 
factor of 15 from the truth on the interval [0,20]. Beyond x — 30, we see what 
appears to be the formation of a spike in the graph. This is best explained by the 
computed values £i (34.858) = 0.00016504... and £i (34.859) = -0.000095463.... 
Since the actual error E 10 is clearly a continuous function, the intermediate value 
theorem implies that Eiq vanishes at some point in the interval [34.858,34.859]. 
On the other hand, the estimated error B\q is clearly an increasing function, so 
it follows that there exist exceptional tiny intervals on which the estimated error 
overestimates the actual error by arbitrarily large factors. However, this kind of 
phenomenon is in general to be expected for Simpson's Rule no matter what sort 
of error estimate is used (apart from one based on the actual error function itself!) 
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Figure 1. Graph of Rio (x) 



4. Understanding the Examples 

Although the reasoning of Section [5] may be easy to follow, the student might 
legitimately wonder how one arrives at useful integral representations such as (|2.2p , 
(|2.6|) , (|2.8p . or (|2. 13[) in the first place. We begin by observing that the examples 
of Section [2] deal with integrands of the form (/(£) — f(0))/t. More generally, if 
/ : [a, b] — ► R has n+1 continuous derivatives, repeated integration by parts applied 
to the equation 

f(t) = f(a)+ [ f'(u)du 

J a 

yields Taylor's formula with the integral form of the remainder: 

/W = E^r/ (j) W+ f { ±^f^\u)du, a<t<b. (4.1) 
— 7 /„ n! 

j=0 J Ja 

If we make the change of variable u = (t — a)s + a and rearrange (|4.1|) . we obtain 



(t - a)™+! 



./("+!) ((i - a ) s + a)ds, (4.2) 
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and hence if / is sufficiently diffcrentiable, Leibniz's rule for differentiating under 
the integral sign yields 

fc /(*)-E"=o(*-«) J 7 0) (a)/i! 



(t - a)" +1 



If a uniform upper bound on the absolute value of the derivative of / of order 
n + k + 1 is known, say 



sup 

a<t<b 



(ri+fe+l) 



(*) 



then we have the upper bound 



— / s k (l-s) n ds 



< M, 



klM 



(n + k+ I)! 

for the absolute value of the left hand side of (|4.3|) . i.e., 



d\ k f(t)-EU(t-*) j f {j) (a>)/j\ 



dt 



(t - a) n 



+ i 



< 



klM 



(n + k + iy: 



(4.4) 



(4.5) 



(4.6) 



It is readily apparent that (|2.5p , (|2.f 0[) . and (|2.I4[) are special cases of (|4.6[) in which 
a = ri = and /(t) = sint, f(t) = — cost, and /(t) = — exp(— t), respectively. 
For an example with n > 0, consider 

I — cos t 



■dt, 



t 2 

which arises in the study of incomplete versions of the Frullani integral (0 p. 470] 
or pH p. 398]) 



cos at — cos pt 
T 2 



dt 



1 — cos f3t 1 — cos at 



f 2 



t 2 



dt 



= (\P\-\a\)-, a,(3eR. 
Repeated application of the product rule for differentiation gives 

d\ 4 l — cost 120(1 — cost) 96sint 36cost 8sint cost 



dt 



t 2 



t 4 



t :i 



t 2 ' 



which again, is troublesome to estimate as it stands. On the other hand, (|4.3[) with 
/(t) = cost, n = 1, a = 0, and k — 4 gives 



whence 



d_ 

dt 

d 
dt 



1 — cos t 

t 2 



I — cos t 



s 4 (I — s) cos(st) (is, 



t 2 



< 



s 4 (I - s)ds 



30' 



with equality at t — 0. More generally, (|4.4p with M = I is satisfied when f(t) 
cost, and so (|4.5[) and (|4.6p give 



I — cos t 



t 2 



< 



s fc (l - s)ds 



(fc+l)(fc + 2) ! 
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with equality at t = when k is even. 



5. Some More Sophisticated Examples 
The inverse tangent integral 

Ti 2 (z)= / -^-dt, 



would appear to be another example amenable to our technique. Lewin [8, chap- 
ter 2] devotes a whole chapter to the study of its properties. The inverse tangent 
integral evidently also attracted the interest of Ramanujan |10j . who among other 
things used it to develop a rapidly convergent series of hyperbolic functions for 
Catalan's constant, Ti 2 (l) ([B p. 807], g]). Formula (|3~2"j) with a = n = and 
/(*) = tan _1 (t) yields 

tan _1 (t) f 1 ds 

(5.1) 



t Jo 1 + sH 2 ' 

After differentiating under the integral and simplifying the resulting expression, one 
arrives at 

d\ 4 tan- 1 (t) f 1 24s 4 (5s 4 t 4 - 10sH 2 + 1) , 

di) — r- = h or^? ds - (5 - 2) 

It is difficult to obtain a good uniform estimate on the size of the integrand in (|5.2p . 
so we return to (|5 . 1 [) . After making the change of variable s = 1/u and employing 
the formula 

poo 

U 

e~ uv cos(vt) dv = — r, u > 0, 

y J u 2 +t 2 

an application of Fubini's theorem [SJ p. 67] gives 
tan _1 (i) f°° 1 udu 



t Ji u u 2 + 1 2 

roo 1 poo 



- wv 







cos(wt) dv du 
du cos(vt) dv 
Ex(w) cos(vt) dv, 

where Ei is the exponential integral (|2.11[) . Thus, 

d\ 4 tan-Vi) f°° a„ , x 
— J 1= v^x{v)cos{vt)dv, 
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and hence by Tonnelli's theorem [6l p. 67], 



tan _1 (t) 



< 



v 4 Ei(u) dv 



■ du dv 



v i e uv dv du 



o 
du 

U 6 



= 4! 
24 

In the same manner, it can be shown more generally that for all nonnegative integers 
k, 

k ■ - . 

v k F,i(v)dv 



tan _1 (t) 



t 



< 



k\ 



V 



with equality at t = when k is even. 

The previous example gives an idea of what is possible when one steps outside 
the framework of Section [5] As another example, suppose one wanted to estimate 
the even derivatives of the tangent function. From [JJ p. 



tant 



sinh 2st 



ds, 



sinh 7rs 

It follows that for all nonnegative integers k, we have 



-it/2 < t < tt/2. 



2A- 

— I tan t = 2 
dt 



As s gets large, we expect 



sinh its 



-tt/2 < t < tt/2. 



(5.3) 



sinli7rs e TS — e _7rs 
to behave like e ( 2 *- 7r ) s for < 2t < tt. In fact, it is easy to prove that the inequality 

sinh 2s\t\ 



< e (2|t|-i0 



sinh its 

holds for — tt < 2t < tt and s > 0. Therefore, by 

2k 



d , 

— ) tanr 
dt 



< 2 



( 2s )2fe e (2|t|-7r)a ds 



we have 

(2k)l 



(TT/2-\t\) 2k + 1 



(5.4) 



for all nonnegative integers k. The estimate (|5 .4[) is remarkably good even when \t\ 
is small, with asymptotic equality in the limit as \t\ — > tt/2—. 

The previous two examples succeeded because we were able to represent the 
desired function as a Fourier (respectively, Laplace) transform of a well-behaved 
function. As a final example, consider the integral 

cb 

t -2K e -ut eK Ein(t) dtj ( 5 5 ) 
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in which < a < b, k > 1, and u > are real parameters, and Ein is the complemen- 
tary exponential integral (|2.12|) . The integral ()5.5j) arises in the solution of certain 
advanced argument difference-differential equations relating to sieves ([2], [3]), and 
as such, it is desirable to be able to compute it accurately for various values of the 
parameters. To study the integrand of (|5.5|) . we let c > and define 

X K (v) := -— I e vz z- 2K e KEin{z) dz. (5.6) 

By Cauchy's theorem, the line integral (|5.6[) is independent of c, and vanishes for 
v < 0. It is not hard to show that A K is continuous and satisfies the delay differential 
equation 

(v^X^v))' =kv- k X k (v-1), v>0, (5.7) 
with boundary condition 

K(v) = ^v K -\ 0<v<l. (5.8) 

r(«) 

Let f{t) denote the integrand of (|5.5p . By Laplace inversion, 



e -vt- 



t X K (v)dv = t- 2K e KEin{t) , t>0, 
and differentiating under the integral sign to obtain 

/W(t) = t- 2K e- ut e KEi ^ = J^iu + v^e-^+^X^dv, 

/(B) ( t ) = ^ t -2« e -«t e « Ei n(t) = _^ °( u + w )5 e -(«+^A K ( w )d w , 

can be justified. But, the delay differential equation (|5.7[) and boundary condi- 
tion (|5.8|) together show that A K is nonnegative, and hence is nonnegative and 
nonincrcasing. It follows that 



sup 

a<t<b 



/W(t) =/ w (a), 



(4)/ 



a considerable simplification. 

6. Conclusion 

Of course, there are other methods for bounding the derivatives of a suitable 
function, the most familiar of which is undoubtedly Cauchy's inequality (see eg. 
[TT1 p. 91]) from the theory of complex variables. One can also consider alter- 
native approaches to error analysis which do not involve estimating higher order 
derivatives. Chapter 4 of Davis and Rabinowitz [5] is devoted to error analysis 
for various approximate integration schemes. In addition to a section on error 
estimates via analytic function theory, alluded to previously, there is a lovely sec- 
tion (see pp. 317-332) describing applications of functional analysis to numerical 
integration and error estimation. 

In this paper, however, we have deliberately focused on the use of integral trans- 
forms arising primarily in the context of real variable theory. The technique can 
be summarized as follows: find a suitable integral representation for the function 
whose derivatives are to be estimated, differentiate under the integral sign, and 
estimate the resulting integral. The technique is hardly new, but is seldom used to 
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its fullest advantage. It is hoped that some of the examples provided herein could 
be used to enrich the discussion of numerical integration in a typical calculus or 
numerical analysis course. 
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